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ABSTRACT 

We present an X-ray spectral and timing model to investigate the broad and variable 
iron line seen in the high flux state of Mrk 335. The model consists of a variable 
X-ray source positioned along the rotation axis of the black hole that illuminates the 
accretion disc producing a back-scattered, ionized reflection spectrum. We compute 
time lags including full dilution effects and perform simultaneous fitting of the 2-10 
keV spectrum and the frequency-dependent time lags of 2.5-4 vs. 4-6.5 keV bands. 
The best-fitting parameters are consistent with a black hole mass of ss 1.3 x 1O 7 M 0 , 
disc inclination of 45° and the photon index of the direct continuum of 2.4. The iron 
abundance is 0.5 and the ionization parameter is 10 3 erg cm s~ 3 at the innermost 
part of the disc and decreases further out. The X-ray source height is very small, 
~ 2 r g . Furthermore, we fit the Fe L lags simultaneously with the 0.3-10 keV spectrum. 
The key parameters are comparable to those previously obtained. We also report the 
differences below 2 keV using the xillver and REFLIONX models which could affect the 
interpretation of the soft excess. While simultaneously fitting spectroscopic and timing 
data can break the degeneracy between the source height and the black hole mass, 
we find that the measurements of the source height and the central mass significantly 
depend on the ionization state of the disc and are possibly model-dependent. 

Key words: accretion, accretion discs — galaxies: active — galaxies: individual: 
Mrk 335 — X-rays: galaxies 


1 INTRODUCTION 


The narrow line Seyfert 1 galaxy Mrk 335 (z = 0.0258) has 
been observed by several X-ray observatories over a range 
of spectral states. In 2006, XMM-Newton found that in a 
high flux state Mrk 335 exhibits very clear double-peaked 
Fe Ka and Fe XXVI Lya emission lines around 6.4 and 7 
keV, respectively ( O’Neill et al.|[2007 Larsson et al.|[2008 |. 
The best-fitting model revealed evidence of relativistically 
blurred inner disc reflection but the narrow peaks at 6.4 and 
7 keV are due to the reflection from more distant material 
such as the molecular torus, and optically thin ionized gas 
Ailing the torus, respectively. In 2007, the source went into 
an extremely low flux state in which the overall flux dropped 
by a factor of « 10 and the spectrum was successfully de¬ 
scribed by either partially covering absorption or blurred re¬ 
flection models ( Grupe, Komossa fe Gallo|2007 Grupe et al. 
20081. In 2009 when the source was in an intermediate flux 


state, the spectral fitting required only the blurred reflection 
from the disc around a rapidly spinning black hole without 


invoking partial covering (Grupe et al. 2012 Gallo et al. 


2013). Recently, Parker et al. (20141 studied the spectrum 


from 2013 observations and constrain the spin parameter of 
the central super-massive black hole to be a « 0.98. 

Alongside time-averaged spectral data analysis, tim¬ 
ing properties of Mrk 335 have been studied. Similar to 
most AGN, it has positive time lags at low frequencies in 
which higher energy bands lag lower energy bands. |Arevalo,| 
McHardy & Summons (2008) found that the sharp cut-off 


in the lags occurred at a frequency close to the bend fre¬ 
quency in the power spectrum which could be described by 


a model of inward propagating fluctuations (e.g., Lyubarskii 


1997 |Kotov, Churazov fc Gilfanov 200T1 |Arevalo fe Uttle~ 


2006). This model assumes the X-ray emission hardens to¬ 


and emissivity profile of Mrk 335 by using the NuSTAR data 


wards the centre, so lower frequency variability in soft emis¬ 
sion at larger radii propagates inwards and modulates the 
higher frequency variability in hard emission at smaller radii, 
thereby producing the positive hard lags. Additionally, many 
AGN with rapid variability show negative soft lags at high 
frequencies, referred to as reverberation lags (i.e. the delay 
between changes in the X-ray continuum and the associated 
echoes from the disc). Since the reflected photons travel an 
additional distance to the disc before being observed, the 
soft reflection variations lag behind the primary hard-band 
variations, with the amplitude of the lag being associated 
with the light travel time between the X-ray source and 
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the disc (see Uttley et al. 2014| for a review). Measuring 
reverberation lags provides a new method to probe the lo¬ 
cation/size of the X-ray source. The first hint of a negative 
lag was seen by McHardy et al. (20071, with the first robust 


detection of the Fe L lags being in 1H0707-495 (Fabian et al. 


2009). After that the Fe L lags have been detected in many 
AGN including Mrk 335 (e.g., De Marco et al.|[2013 l. Fur¬ 
thermore, Kara et al. (2013) found that the Fe K lags are 


also present in Mrk 335 during the high flux state. These 
lags disappeared once the source went into an intermediate 
flux state. Combining the high and medium flux state data 
showed lags are comparable to those from the high flux state 
alone, demonstrating that the lag signal is dominated by the 
high flux data. 

The interplay between the disc reflection and propagat¬ 
ing fluctuations seems to provide a complete explanation of 
time lags in AGN. This framework is preferred over alterna¬ 
tive models (e.g., reflection from distant materials proposed 
by Miller et al. 2010) because it can explain, for example, 


the absence of energy dependent Fe K lags at low frequencies 
(Kara et al. 2013) as the low- and high-frequency lags are 


generated by two distinct mechanisms. The reverberation 
occurs at shorter timescales compared to the inwards prop¬ 
agating fluctuations so they dominate at higher frequencies 
while the propagation lags take over at lower frequencies. 

In this paper we analyse the high flux data from XMM- 
Newton during observations in 2006. We focus on self- 
consistently modelling time-averaged and and lag-frequency 
spectra of Mrk 335 in the lamp-post geometry scenario (i.e., 
the disc is illuminated by an axial isotropic point source). We 
first investigate the lags between the hard bands, 2.5-4 vs. 4- 
6.5 keV, before shifting the lags to two softer bands, 0.3-0.8 
vs. 1-4 keV. The photon trajectories are traced along null 
geodesics in the Kerr metric between the source, the disc and 


the observer i 

Bardeen, Press & Teukolsky 

1972 

Cunning- 

ham 1975| |Karas, Vokrouhlicky & Polnarev 

1992 

Reynolds 


We assume the disc to be geometrically thin, optically thick 


(Shakura & Sunyaev 1973 I and apply the xillver model 
( Garcia fe Kallman|2010 Garcia et al.|2013 1 to deal with the 
X-ray reprocessing by the disc via fluorescence, Compton¬ 
scattering and photoelectric absorption. The ionization state 
of the disc is determined given the disc density and the 
illuminating flux. The total spectrum is divided into two 
components: Power-Law dominated Component (PLC) and 
Reflection Dominated Component (RDC). In principle the 
RDC lags behind the PLC since it is dominated by reflected 
photons which take a longer time before being observed. The 
reflection responses to the primary variations (so-called re¬ 
sponse functions) are used to compute their delays in the 
Fourier frequency domain, indicative of reverberation lags. 
We parametrise the positive continuum lags by a power-law, 


Emmanoulopoulos et al. (2014). 


Nevertheless, each energy band contains both power- 
law and reflection components. The measured lags then are 
diluted as there is some continuum emission in the RDC 


and some reflection emission in the PLC (Wilkins & Fabian 


2013). Recently Cackett et al. (20141 investigated the effects 


of key parameters such as the source height, inclination angle 
and the black hole mass on time lags. They applied the re¬ 
flected response fraction to the RDC to dilute the Fe Ka lags 
of NGC 4151, and assume no contribution of reflection to the 


PLC. Emmanoulopoulos et al. (20141 successfully fitted the 


Fe L lags of 12 AGN by computational modelling the re¬ 
verberation lags and employing a power-law for the positive 
lags. They, however, assumed the contamination of the con¬ 
tinuum in the RDC and no reflection flux in the PLC. Here 
we include the full contamination between cross-components 
(soft in the hard band and hard in the soft band) so that 
full dilution effects are taken into account and, as a result, 
the lags are further suppressed. As our model predicts time- 
averaged energy spectra and frequency-dependent reverber¬ 
ation lags, the low frequency positive lags are modelled sep¬ 
arately by a phenomenologically motivated power-law. All 
fitting was performed using ISIS ( Houck fe Denicola|2000 ). 

The remainder of this paper is organised as follows. We 
describe the XMM-Newton observations and data reduction 
in Section [2] The reverberation model based on ray trac¬ 
ing simulation is presented in Section [3] The model is in¬ 
vestigated for cases in which there is a density gradient in 
the ionized disc. In Section [4] we present the model pa¬ 
rameters for Mrk 335 and perform simultaneous fitting of 
the time-averaged and lag-frequency spectra. The conclu¬ 
sions are presented in Section [5] We discuss the differences 
we find when using the reflionx model (George & Fabian 


|1991[ |Ross, Fabian fc^Yo ung 1999s |Ross fc “F abian 2005|) in 

the Appendix. 


2 OBSERVATIONS 


Mrk 335 was observed by XMM-Newton (Jansen et a l. 
20011 on 2006-01-03 for 133 ksec. The time lag and time- 


averaged spectra were produced in a standard way, following 
the data reduction of Ka ra et al.| |2013| using the XMM- 
Newton Science Analysis System version 14.0.0. The data 
were screened for background flares, and filtered using the 
conditions PATTERN <= 4 and FLAG == 0, leaving 120 ksec 
of good data. The source spectrum and light curves were 
extracted from a circular aperture of radius 35”, and back¬ 
ground counts were extracted from a similarly sized but off¬ 
set aperture. Background subtracted lightcurves were pro¬ 
duced using the lccorr tool with 10 s time bins. The time 
lags were computed from the phase of the average cross 
power spectrum in the standard way, as were the error bars 


on time lags, following Nowak et al. (19991. The time lags 


that we obtain are consistent with those previously pub¬ 
lished in the literature. 


3 SPECTRAL AND TIMING MODEL 

This section summarises essential equations and theoretical 
framework relating to our physical reverberation model. We 
will first describe the ray tracing simulation in Section 3.1 
followed by the response function which is used to produce 
the full-contaminated light curve in Section 3.2. In Section 
3.3 we investigate the reverberation lags and discuss the 
effects of dilution, ionization, disc density, and disc outer 
radius on the time lags and energy spectrum. 
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3.1 Ray tracing and X-ray reprocessing 


The Kerr metric describing the spacetime around a ro¬ 
tating black hole is characterised by the black hole mass, 
M, and the angular momentum, J, parametrised in the 
form of the specific angular momentum or spin parameter 
a = J/M. Throughout the ray tracing simulation we use 
Boyer-Lindquist coordinates and measure the distance and 
time in gravitational units, r g = GM/c 2 and t g = GM/c 3 , 
respectively. The disc is assumed to be a standard geometri¬ 
cally thin, optically thick disc (Shakura & Sunyaev 1973) 


extending from the innermost stable circular orbit, r ms , 
to 400r g . We assume the observer is at 1000r 9 from the 
black hole and consider a lamp-post geometry in which an 
isotropic point source is positioned at height h along the 
symmetry axis. 

To compute the direct and reflection spectra, an adap¬ 
tive step-size method is used to trace all photon trajectories 
along the null geodesics between the source, the disc and the 
observer. The source spectrum is a powerlaw with flux den¬ 
sity given by F S (E S ) oc Ej a , where E s is the photon energy 
and a is the energy index. All conserved quantities are the 
total energy, E, the angular momentum, L, and the Carter 


1975 


1997 


constant, Q (Carter 19681. The equations of motion of a 
photon are (Bardeen, Press fc Teukolsky|1972 [Cunningham 


Karas, Vokrouhlicky & Polnarev)|1992[ [Fanton et al, 

Reynolds et al.||1999 1 


where 


dr 

dA 

= ±(V r ) 1/2 , 



(1) 

dd 

dA 

= ±(Ve) 1/2 , 



(2) 

^d\ 

= —(aE — LI sin 2 

e) + 

aT/A, 

(3) 

dt 

dA 

= —a(aE sin 2 9 — 

L) + 

(r 2 + a 2 )T/ A, 

(4) 

E = 

2 i 2 2/i 

r + a cos 0, 



(5) 

A = 

r 2 + a 2 — 2 Mr, 



(6) 

T = 

E(r 2 + a 2 ) — La, 



(7) 

K = 

T 2 - A[/xV + (I 

/ — aE) + Q], 

(8) 

Ve = 

Q — cos 2 9[a 2 (g 2 - 

- E 2 ] 

) + L 2 / sin 2 0\, 

(9) 


Note that A = r//r where r and g are the proper time and 
the rest mass of the test particle, respectively. 

We perform photon path integrals numerically in paral¬ 
lel on NVIDIA K20 Graphic Processing Unit (GPU) cards 
on the BlueCrystal Supercomputer at the University of Bris¬ 
tol. Those photons from the source that arrive at the ob¬ 
server directly are given the direct spectrum, 

Fn(E 0 ) = gZ +1 E~ a dE 30 /dS 0 , (10) 


where dE so is the solid angle corresponding to the surface 
area dS 0 on the observer sky and g so is redshift of the source 
in the observer’s frame. The redshift, g so is dehned as 


Qso — 


Vs 


Po • Uo 
Ps U s ’ 


( 11 ) 


where p 0 ( s ) is the 4-momentum of the observed (source) 
photon and u D ( s ) is the 4-velocity of the observer (source). 
Furthermore, we divide the disc into small drd(f>- bins. 


X-ray Reverberation in Mrk 335 


The incident flux per unit area of each bin is then given by 

Fd{E d ) = g^E^dZsj/dSa, (12) 

where g s d is redshift between the source and the disc, dSd is 
the unit surface area in the disc frame and E s d is the solid an¬ 


gle corresponding to the disc element drdcj) (see Ruszkowski 


2000 |Dovciak, Karas fc Yaqoob 2004; Chainakun & Young 


2012 for a full expression). The corresponding ionization 
parameter is defined as 


£(r, v?) = 


4nFt(r, ip) 
n(r) ’ 


(13) 


where Ft (r, i p) is the total illuminating flux in that bin. The 
disc density is assumed to follow a power-law form: n(r) oc 
r~ p . We apply the xillver model (Garcia fc Kallman|2010 


Garcfa et al. 2013I to obtain the reflection spectrum allowing 


the iron abundance, Aj? e , and the photon index, F, to be 
model parameters. A comparison with the reflionx model 


( George fc Fabian||1991 Ross, Fabian & Young 1999; Ross 


& Fabian 20051 will be given in the Appendix. Note that 


the photon index relates to the energy index via F = a + 1. 

To calculate the illumination of the disc we trace ~ 
5 x 10 7 photons from an isotropic point source to the disc 


as described by Chainakun & Young (20121. To trace the 


photons from the disc to the observer, we integrate photon 
paths backwards through 8192 x 8192 pixels on the observer’s 
sky to the disc (Fanton et al. 19971. The reflection spectrum 
is blurred once it is transferred to the observer’s frame. The 
total time-averaged spectrum is the sum of the direct and 
blurred reflection components. 


3.2 Response function and full, contaminated 
light curve 

The response function, ip(t), which is used to produce time 
lags is the count rates of reflection photons as a function 
of time after the primary X-ray variations. However the re¬ 
flection spectra from the disc could be different, resulting in 
different response profiles for the soft and hard energy bands, 
and ipH(t), respectively. Examples are shown in Fig. |T] 
for the responses of the 0.3-1 keV and 5-7 keV bands to an 
instantaneous flash of a point source at h = 5r g . Other pa¬ 
rameters are iron abundance Ape = 1, photon index T = 2, 
spin a = 0.998, inclination i = 30° and ionization state at 
the innermost stable circular orbit £ ms = 10 3 erg cm s _1 . 
For a disc with ionization gradient the relative changes in 
spectra between each energy band are different correspond¬ 
ing to various ionization states on the disc. The differences in 
the responses are the result of the different energy bands se¬ 
lected. This implies the choice of the energy bands can affect 
the response function and hence the time lags. Throughout 
this paper we normalise the areas of the response functions 
to 1 and apply soft and hard reflected response fractions, 
Rs and Rh , respectively, to determine the relative levels 
of their responses. These fractions also determine the levels 
of contamination in both bands which cause dilution effects 
and, as a result, the measured lags is less than the intrinsic 
lags. 

Assuming a(t) is the primary X-ray continuum flux as 
a function of time, we derive the full, contaminated soft and 
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Time (tg) 

Figure 1. Response functions of the 0.3-1 keV (black, solid line) 
and 5—7 keV (red, dashed line) bands normalised so that their 
underlying areas equal 1. The relative changes of spectra between 
both bands for a disc with ionization gradient are not only in 
amplitude but also in shape which could lead to differences in the 
profiles. See text for more details. 


hard band light curves, S(t) and H(t ), respectively, as 

S(t) = as(t) + Rs [ — t')dt' (14) 

Jo 

and 

H(t) = a,H{t) + Rh f a(t')ipH{t — t')dt'. (15) 

J o 

The source flux can be varied to produce several ioniza¬ 
tion states on the disc. However, for simplicity, we select 
only one ionization profile each time when we do Fourier 
transform of the response function. In other words we as¬ 
sume the ionization gradient does not change over a pe¬ 
riod of time t for the primary variation a(t). Since the di¬ 
rect continuum follows the primary variations, as(t) and 
an(t) are interpreted as the contribution from the direct 
power-law components in soft and hard bands, respectively. 
The integral terms represent the response of the disc to 
earlier variations that contribute to the reflected flux. Rs 
and Rh are, respectively, the soft and hard reflected re¬ 
sponse fractions defined, in the same way as other authors, 
as (reflection flux)/(continuum flux) indicative the contam¬ 
ination ratio between cross-components. The fraction Rs 
is one of our model parameters while Rh is then measured 
from the corresponding spectra. Equations (14) and (15) are 
normalised such that Rs(h) = 1 corresponds to equal con¬ 
tribution of the direct and reflection components in the soft 
(hard) band. 


3.3 Frequency-dependent time lags 

The nature of reverberation lags in AGN has been investi¬ 


gated through a wide range of parameters (Wilkins & Fabian 


2013 Cackett et al. 2014 Emmanoulopoulos et al. 2014). 


These lags are sensitive to the source’s height as they relate 
to the light-crossing timescales between the source and the 
disc. If the black hole mass increases, the lags are system¬ 
atically shifted to lower frequencies and higher amplitudes 


(De Marco et al. 20131. The spin parameter and the incli¬ 
nation angle, however, seem to have less effect on the lags. 
See Cackett et al. (2014) for a detailed comparison of lags 
for different geometries of the accretion flows. Throughout 
this section we select the soft and hard bands to be 0.3- 
0.8 keV (RDC) and 1-4 keV (PLC) bands, respectively. 
When the power-law flux contaminates the RDC, the re¬ 
flected response fraction in the RDC decreases and time lags 
decrease. On the other hand, when the reflection flux con¬ 
taminates the PLC, the reflected response fraction in the 
PLC increases and time lags decrease. The more contam¬ 
ination between cross-components the more the time lags 
are diluted. Current literature such as Cackett et al. (20141 
and Emmanoulopoulos et al. (20141 include the direct flux 
in the RDC but assume no reflected flux in the PLC. Here 
we produce the full, diluted time lags by also including the 
reflection component in the PLC. 

We investigate further the case of a density, and hence 
ionization, gradient in the disc and the dilution effects of 
contamination between the hard and soft bands. We follow 


Nowak et al. (1999) to calculate the frequency-dependent 


time lags by taking the Fourier transform of the full, con¬ 
taminated light curves, S(t) and H(t). Their cross-spectrum 
is 


C(f) = S*(f)H(f). 


(16) 


The argument of the cross-spectrum gives the phase differ¬ 
ence, <(>(/) = argC(/), which relates to time lags via 


>-(/) = 


Hf) 

27 Tf ■ 


(17) 


In our definition, a negative lag means the soft band lags 
behind the hard band. The phase difference is defined over 
—7r to 7r, so phase wrapping occurs at / ^ l/2r where the 
profile fluctuates around 0. 

The necessity of including the full components in the 
soft and hard band light curves is emphasised by [Wilkins 
& Fabian (2013). This is achievable by our model. The 


variations in flux-contamination can be physically produced 
when, for example, the ionization state of the disc changes 
in response to changes in the source luminosity. Fig. [^rep¬ 
resents the reflected response fraction and lag spectra for 
the source at 5 r g whose luminosity is varied to produce low, 
medium, and highly ionized discs. Fixing Rs = 1.2, the mea¬ 
sured Rh is 0.17, 0.53 and 0.9, respectively. Other param¬ 
eters are i = 30°, Ap e = 1 and T = 2. Even though the 
light-crossing distances are constant, the time lags decrease 
with increasing ionization state, and for the highly ionized 
disc the lags are very small. This is because when the disc 
becomes more ionized, the reflection spectra become flatter 
so Rh increases and, as a result, the lag amplitude decreases. 

According to the model, the ionization state depends 
not only on the total incident flux, but also on the disc 
density which we assume to vary radially as n(r) oc r~ p . 
Fig. [3] shows how the ionization state and lags depend on the 
density index, p. For comparison, we fix the ionization state 
at the innermost radius, £ ms , to be 10 5 erg cm s -1 . As the 
inner disc is more strongly illuminated than the outer parts, 
the ionization state should decrease further out. Decreasing 
p increases the slope of ionization profile covering broader 
ranges down to very low states at the outer radii. This leads 
to the smaller Rh and larger reverberation lags. 
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Figure 2. Top panel: The reflected response fraction for the cold 
(black line), medium-ionized (red, dashed line) and highly-ionized 
(blue, dotted line) disc. We allow the ionization states, given a 
value at the inner radius, to vary across the disc. The measured 
Rh for fixing the Rg = 1.2 are 0.17, 0.53 and 0.9 respectively. 
The soft and hard bands are highlighted by yellow and blue back¬ 
ground shading, respectively. Bottom panel: The corresponding 
time lags as a function of Fourier frequency. As the disc becomes 
more highly ionized the amplitude of the time lag decreases. 

Although elsewhere in this paper we assume the outer 
radius of the disc r 0 = 400r g , here we investigate the effects 
of relaxing this assumption. Time lags and the reflection 
spectra for variable r 0 are presented in Fig. [4] When the disc 
is extended further out, there are more regions for reflection 
associated to longer time delays. Therefore the averaged re¬ 
sponse time increases and hence time lags increase. We can 
see the maximum lag difference between r 0 = 20 and 400r g 
is « 5 t g . However, their reflection spectra are almost identi¬ 
cal which agree with the theoretical framework in which the 
prominent fluorescent lines (e.g., the broad Fe Ka line) are 
produced at the inner disc so they are almost independent 
on r 0 . 

It is clear that while the parameter r 0 has less effects on 
the averaged spectrum, it does affect time lags. There is also 
the possibility that the response fractions in soft and hard 
bands increase with r 0 because of an extended reflecting 
area. Since the outer regions have a lower ionization state 
than the inner regions, the ratio Rs/Rh is larger for the 
outer part of the disc. If this is the case, the maximum lag 
differences are larger than those in Fig. 5, which assumes 
equal Rs = 1 for all cases. However, to limit the number of 
parameters we restrict our consideration to r 0 = 400r g . We 



Frequency (1/tg) 

Figure 3. Top panel: ionization state of the disc for different den¬ 
sity profiles, n(r ) oc r~ p . Bottom panel: The corresponding time 
lags. In the scenario simulated here, discs with stronger density 
gradients, i.e., larger p, show shorter lags because the disc can 
remain highly ionized out to larger radii. 

note that r 0 could affects time lags, and could be included 
as a variable parameter in the future. 

4 FITTING MRK 335 
4.1 Model parameters 

We consider a maximally spinning black hole and assume the 
accretion disc extends from r ms to 400r g . The cut-off energy 
of the X-ray continuum is fixed at 300 keV. Our model (so- 
called revb model) is coupled to the xillver model that has 
variable angles of reflection X-rays with respect to the disc’s 
normal between 5° — 85°. For simplicity only one value of 
inclination « 50° is selected but this could be relaxed in the 
future. We produce the grid of parameters shown in Table [T] 
which consists of 40,320 grid cells covering possible values 
that have been found by previous studies. Each grid cell cor¬ 
responds to a given set of parameters: the source height ( h ), 
disc inclination angle (i), photon index (r), iron abundance 
(Ai?e), ionization state at the innermost stable circular or¬ 
bit Rms), the disc density index ( p ) and the soft reflected 
response fraction (Rs)- A model without the response frac¬ 
tion is very challenging and depends upon many factors (e.g. 
the efficiency of reflection or the normalisation of the mod¬ 
els). Although ray tracing simulation allows us to count the 
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Figure 4. Outer radius variation for the low ionized disc (£ < 
10 erg cm s —1 and A = 1) illuminated by an on-axis point source 
with h = 5 r g and T = 2. Other parameters are a = 0.998, 9 = 30° 
and Rg = 1. Top panel: The reverberation lags for five values 
of the outer radius. Bottom panel: the corresponding reflection 
spectra. While changing the outer radius does not significantly 
affect the shape of the reflection spectrum, discs with a large 
outer radius show longer time lags. 


number of hard photons that escape to the observer and that 
irradiate the disc, the proportionality coefficient of the re¬ 
flection and the incident flux cannot be determined without 
further assumption of such as the albedo of the reflecting 
disc. This makes calculating the exact value of the response 
much more complicated. For this reason we allow the frac¬ 
tion Rs to vary and measure the corresponding Rh . In other 
words the Rs is one of the model parameters which once se¬ 
lected will determine a unique value of Rh- Note that our 
model grids extend over a large, multi-dimensional param¬ 
eter space which is computationally difficult to interpolate 
over, and when fitting data will result in a complex \ 2 space 
that will have many local minima. We therefore have cho¬ 
sen to step through the entire grid of models and allow, for 
example, the black hole mass and normalisation factors to 
be free parameters. 

The models that we have computed predict both the 
time-averaged spectrum and reverberation lags, so they can 
be simultaneously fit to the data. Since our reverberation 
lag model does not explicitly compute the positive lags, 
we employ a power-law (powerlaw model) in the form 
of t(/) = N(cf)~ s to produce the positive lags where 
N 6 [0,10 6 ], s £ [0,6] and c = 10 4 is a scale factor 


Parameters 

Value 

h ( r g ) 

2, 3, 4, 5, 6, 7, 8 

i (°) 

15, 30, 45, 60 

r 

1.8, 2, 2.2, 2.4 

A-Fe 

0.5, 1, 2 

log £ms( erg cm s -1 ) 

1, 2, 3, 4, 5 

P 

0, 1, 2 

Rs 

0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 


Table 1. Grid of parameters for which the REVB model has been 
computed. 


of the frequency /. Therefore the total timing model is 
revb + powerlaw. The time in gravitational units is trans¬ 
formed to physical units for the variable black hole mass 
between 10 6 — 1O 9 M0. For the spectral model we add two 
narrow Gaussain lines at « 6.4 and 7 keV responsible for the 
refection from very distant material as discussed in|Q’NeiTf| 


et al. (2007) and Larsson et al. (20081, and necessary ad¬ 


ditional components when modelling the broad iron line. 
These narrow features do not affect the lags as they do not 
vary on the timescales of the inner disc reflection. All spec¬ 
tral components are modified by galactic absorption with 
Nh = 3.99 x 10 2l) cm -2 . The total spectral model is pro¬ 
duced via TBABS ® (revb + GAUSS 1 + GAUSS2). 


4.2 Results 


We begin by fitting the 2-10 keV spectrum and the 2.5-4 
vs. 4-6.5 keV time lags in high flux state of Mrk 335. The 


fitting procedure is performed in ISIS (Houck & Denicola 
2000) using the SUBPLEX minimisation method. In ISIS we 
load the time-averaged XMM-Newton EPIC-pn spectrum, 
background and responses as the first data set. We then de¬ 
fine a second data set which represents the time lag vs. fre¬ 
quency, rather than the standard counts vs. wavelength. For 
a given set of parameters our reverberation model produces 
the time-averaged spectrum which is fit to the first data set 
and the time lag spectrum that is, simultaneously, fit to the 
second data set. The reported \ 2 values of fit statistics are 
for both data sets combined. We have been able to achieve 
excellent fits without having to change the relative weight¬ 
ing of the spectroscopic and timing data sets, although this 
can be done if required. 

The best fit spectral and timing profiles are shown 
in Fig. [5] The best-fitting parameters are listed in Ta¬ 
ble [2] The model provides excellent fits to the data with 
X 2 /d.o.f. = 1.026. The X-ray source is at 2r g on the sym¬ 
metry axis. The time lags are well-fitted for the black hole 
mass M « 1.3 x 10 7 Mq. The accretion disc has an inner 
ionization parameter £ ms = 10 3 erg cm s _1 , A_p e = 0.5 and 
p = 0. The inclination angle is 45° and the soft reflected 
response fraction is 0.3. We find the photon index of the 
continuum X-rays is 2.4 steeper than those of normal AGN. 
However, the steep photon index has been found before in 
Mrk 335 by such as |Gondoin et al. (|2002[) and Wilkins & 


Gallo (2015). The fits are poorer without two narrow Gaus- 


sians at « 6.4 keV (Fe Ka) and 7 keV (Fe XXVI). While the 
first line can be interpreted as the result of reflection from 
distant, cold material such as a molecular torus, the sec- 
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Energy (keV) 



Figure 5. Data and residuals (blue points) from simultaneously 
fitting the REVB xillver to the high flux state spectra (top 
panel) and time lags between 2.5—4 and 4-6.5 keV bands (bottom 
panel). The overall model is shown in red, which is the sum of 
individual model components (shown in green). The individual 
spectral components are the continuum, blurred reflection, and 
two iron lines (from top to bottom, respectively) and the individ¬ 
ual lag components are the power-law and reverberation model 
(from top to bottom, respectively). Note that there is clearly a de¬ 
generacy between the the power-law, which may be attributed to 
the propagation of fluctuations in the disc, and the reverberation 
signature. 


ond line may originate in hot, optically thin material such 
as the hot gas filling that torus. Note that the analysis of 


O’Neill et al. (20071 required these two additional narrow 


line components in addition to the relativistically broadened 
reflection component. These narrow features from distant 
reflection do not vary on the timescales of the inner disc re¬ 
flection so including them in the spectral model do not have 
any effect on the lags. Therefore the model can provide a 
self-consistent fit to both spectrum and time lags within the 
same energy band. 

Having successfully fitted a physical reverberation 
model to the 2-10 keV spectrum and time lags of 2.5-4 vs. 
4-6.5 keV bands, we now extend to fit the spectrum down to 
0.3 keV with the Fe-L lags of 0.3-0.8 keV vs. 1-4 keV bands. 
To deal with the soft excess at ss 0.3-2 keV which is often 
seen in AGN, we replace two Gaussians with more physi¬ 
cal models, unblurred xillver and vmekal, responsible for 
distant reflection from cold medium and ionized diffused- 
gas, respectively. The warm absorber as reported by some 


Component 

Parameter 

Value 

TBABS 3 

Ahj(xl0 20 cm- 2 ) 

4.0 f 

REVB 3+t 

H r g) 

2 


i^) 

45 


r 

2.4 


Ap e 

0.5 


log ( erg cm s“ 

- 1 ) 3 


P 

0 


Rs 

0.3 


Norml 

7.6 x 10 5 


log M(M 0 ) 

7.1 

POWERLAW 4 

s 

0.82 


Norm2 

151.8 

GAUSS I s 

Area (x 10 -6 ) 

9.81-f 


Centre (keV) 

6.4f 

GAUSS2 3 

Area (x10 -6 ) 

4.80-f 


Centre (keV) 

7.0 f 

X 2 /d.o.f. 


1.026 


Table 2. Model parameters for fits to the 2—4 keV spectrum and 
time lags between 2.5—4 and 4—6.5 keV bands of Mrk 335 dur¬ 
ing the high flux state. The model components, their parameters 
and the parameter values are listed in Columns 1, 2 and 3, re¬ 
spectively. The superscript s refers to spectral parameters while 
t refers to timing parameters. The superscript / identifies the pa¬ 
rameters which are fixed. Line energies refer to the rest-frame. 


previous studies (e.g. Longinotti et al. 2013I is modelled 
by XSTAR ( Kallman fe Bautista||200T I . Finally, a blackbody 
component with a maximum temperature fixed at 0.05 keV 
is added and is corresponding to the disc radiation which is 
constant over time. Therefore, all components we include so 
far should have no (or very little) effect on the reverberation 
lags. 

The model fitting to 0.3-10 keV spectrum and Fe L 
lags still provides good fits to the data with x 2 /d.o.f. = 1.15 
(see Fig. |6|. Previous studies on the RGS spectrum (e.g., 
Grupe et al.|2008 Longinotti et al.|2008 \ found the atomic 


emission lines from the ionized plasma play an important 
role in the soft bands. In the vmekal model, we allow the 
abundances of O, Si, S and Fe atoms to vary between 0.2- 
5, and tied the abundances of other atoms together. The 
important parameters are listed in Table[3] Most of the revb 
parameters are comparable to what previously obtained by 
constraining the lags in harder energy band except £ ms = 
100 erg cm s -1 , p = 2 and Rs = 0.4. The difference of Rs 
can be understood as the selected soft bands are different. 
The reverberation lags in this case are negative since the soft 
band is more dominated by the reflection flux comparing to 
the hard band. The lags scale with the mass M « 2 x 10 7 Mq . 
The powerlaw components seen in both cases are quite 
different. This is because the continuum lags are dependent 
on the selected energy band as same as the reverberation 
lags. 


4.3 Discussion 

Simultaneously fitting the 2- 10 keV spectrum and time lags 
between 2.5-4 keV and 4-6.5 keV is straightforward for two 
main reasons. Firstly, the lags are in the same energy band 
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Figure 6. Data and residuals from simultaneously fitting the 
REVB (g> xillver to the 0.3—10 keV spectra (top panel) and the 
time lags between 0.3-0.8 vs. 1-4 keV band representing the 
Fe L lags (bottom panel). The individual model components are 
shown in green. In the top panel, the spectral components, folded 
through the instrument response, are the continuum, blurred re¬ 
flection due to reverberation, and xillver neutral reflection (from 
top to bottom, respectively, at high energies), a thermal compo¬ 
nent that exceeds the blurred reflection component around 1 keV, 
and a blackbody component that only contributes significant flux 
< 0.5 keV. In the bottom panel, the lag components represent 
power-law and reverberation (from top to bottom, respectively). 
In contrast to the case shown in Fig. [5] the power-law and re¬ 
verberation lag components have opposite signs so are easier to 
separate. 


as the time-integrated energy spectrum. Secondly, two addi¬ 
tional narrow features « 6.4 keV and 7 keV from the distant 
reflection do not affect the lags because they do not vary on 
the reverberation timescales. However, shifting the lags to fit 
the Fe L band is much more difficult due to the soft excess 
at « 0.3-2 keV and the choices of appropriate additional 
components/models. One important reason why the revb 
model alone does not fit the soft band well is that while 
the iron abundance was allowed to vary the abundances of 
other elements were set to the Solar value. Allowing vari¬ 
able abundances of other atoms in the revb model should 
significantly improve the fits in the 0.3-2 keV band but this 
will make changes in the soft response fraction and rever¬ 
beration lags. Unfortunately these changes are beyond the 
scope of our present paper, although the change to the soft 
band spectrum could be significant (cf. the discussion in 


Component 

Parameter 

Value 

TBABS 3 

ATrfxlO 20 cm -2 ) 

4.0^ 

XSTAR s 

Norm 

1.7 x 10 21 


log 5ms ( erg cm s _1 ) 

1.85 

REVB s+t 

h(r g ) 

2 


i (°) 

45 


r 

2.4 


Af £ 

0.5 


log 5ms ( erg cm s _1 ) 

2 


P 

2 


Rs 

0.4 


Norml 

5.6 x 10 s 


log M(M q ) 

7.3 

POWERLAW 4 

S 

2.92 


Norm2 

419.6 

BLACKBODY s 

T (keV) 

0.05 


Norm 

4.0 x 10" 4 

unblurred xillver s log £ ms ( erg cm s 1 ) 

of 


Afb 

0.7 

VMEKAL S 

T (keV) 

0.51 


Ao 

2.8 


Asi 

0.2 


As 

0.2 


Afb 

0.9 


-Mother 

1.8 

X 2 /d.o.f. 


1.15 


Table 3. Model parameters for fits to the 0.3-10 keV spectrum 
and Fe L lags between 0.3-0.8 vs. 1-4 keV bands of Mrk 335 dur¬ 
ing the high flux state. The model components, their parameters 
and the parameter values are listed in Columns 1, 2 and 3, respec¬ 
tively. The superscript s, t and f refers to spectral parameters, 
timing parameters and fixed parameters, respectively. 


the Appendix on differences between reflection model with 
the same abundances). Instead we model the distant reflec¬ 
tion using (unblurred) xillver, and the soft excess with a 
VMEKAL plus blackbody component at the maximum tem¬ 
perature 0.05 keV. These additional components do not in 
principle significantly affect time lags. Future high resolution 
spectroscopy with Astro-H and Athena will be particularly 
helpful in disentangling the spectral complexity of the soft 
X-ray band. 

We find that one warm absorber is needed for the good 
fit below 2 keV. The unblurred xillver contributes mainly 
to the narrow line at « 6.4 keV but very little to the rest of 
spectrum, so using a Gaussian is still a good approximation. 
The plasma temperature obtained with the vmekal model 
is low (0.51 keV) so it does not provide a perfect fits to the 
« 7 keV weak Fe XXVI line. However, it does remarkably 
improve the fits below 2 keV. The requirement of xillver 
(or « 6.4 keV Gaussian line) and vmekal suits the frame¬ 
work of the reflection from distant cold-molecular torus and 
the ionized gas filling the torus dO’Neill et al. 2007). This 
is a plausible model to fit the 0.3-10 keV spectrum without 
invoking the components that might significantly affect time 
lags. 


The soft excess can also be modelled by either a blurred 
reflection or a blackbody component with a temperature « 
0.1-0.2 keV ( Gierlinski fc Done|2004 Crummy et al.|2006 |. 
Although this temperature is argued to be too hot for the 
thermal emission from a standard disc, it is possible that 
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the inner disc under strong light bending effects is heated 
by the returning radiation, causing the disc to emit the soft 
X-rays. Grupe et al. (2008) successfully modelled the high 
flux state of Mrk 335 using two power-laws, either two or one 
corresponding blurred inner-disc reflection, plasma emission 
by mekal model and some further emission lines. Recently 


Wilkins & Gallo (2015) fitted the same spectrum profile but 


between 0.5-10 keV using the blurred and unblurred com¬ 
ponents modified by the warm absorbers. Our blackbody 
component is negligible if we fit the 0.5-10 keV spectrum 
and time lags {x 2 /d.o.f. « 1.1) because it contributes to the 
flux < 5 keV. Using MEKAL rather than vmekal gives ac¬ 
ceptable but poorer fits (\ 2 /d.o.f. « 1.3). Keeping in mind 
the complexity of the soft excess (e.g. its uncertain origin 
and variation-timescales), our current model is more self- 
consistent if the spectrum and the lags are fitted above the 
2 keV band. 

Moreover, we find that the model still provides a good 
description of the data (x 2 /d.o.f. < 1.2) for all disc density 
indices (p = 0, 1 and 2), but the source height and the black 
hole mass can be different. Of the values we investigated, 
the source is at « 2 — 3r s on the symmetry axis and the 
acceptable black hole mass is « 1.3 — 2.5 x 10‘ Mg depend¬ 
ing on a given density profile of the disc. Our constrained 
mass agrees with the estimated mass that follows the mass¬ 
scaling laws (e.g. De Marco et al.|2013 |. It is in good agree¬ 
ment with the masses reported by Grier et al. (20121 and 


Emmanoulopoulos et al. (2014) derived from optical and X- 


ray reverberation, respectively. The source height is slightly 
lower than Emmanoulopoulos et al. (2014) where h = 3.1r 9 
was found by fitting only the Fe-L lags and neglecting the 
contamination in the hard band. Such low height means the 
corona is very compact even in the high flux state of this 
source. It would be interesting to see how the model param¬ 
eters change when the source is in a low flux state. 

Simultaneously fitting spectral and timing data may 
help break degeneracies between the continuum source 
height and black hole mass. The source height determines 
the illumination pattern and, given the disc density profile, 
the corresponding reflection spectra. The black hole mass 
sets the light-crossing timescale of a gravitational radius, 
and the source height gives the associated reverberation 
lags. Only by combining spectroscopy and timing products 
are both the black hole mass and source height constrained. 
However, different disc density profiles may result in differ¬ 
ent estimates of the source height and black hole mass. This 
suggests that a better understanding of the radial density 
distribution in the disc is needed which in turn allows a more 
accurate estimation of the ionization gradient corresponding 
to a particular illuminating X-ray source. Otherwise there 
will be some degeneracy between the source height and cen¬ 
tral mass produced by different density profiles that still 
provide excellent fits to the data. Fitting energy-dependent 
time lags (e.g., the full lag-energy spectrum) might break 
these degeneracies, which will be investigated in future pa¬ 
pers. 

Finally, it is important to note that even considering 
reverberation between harder X-ray bands is satisfactory 
without taking into account the uncertainty of the “soft ex¬ 
cess” below 2 keV (or, crucially, how it varies), the sign 
of the time lag caused by putative fluctuations propagat¬ 
ing through the disc and reverberation due to light travel 


time are the same (see the green lines in the bottom panel 
of Fig. [5j| . This possible degeneracy will be addressed in a 
future paper by self-consistently modelling the propagating 
fluctuations. This also points towards joint spectral and tim¬ 
ing modelling being more straightforward, from a theoretical 
perspective, in the Fe K band which is observationally more 
challenging but ideally suited to the Athena X-ray observa¬ 
tory (e.g., Dovciak et al.|2013). 


5 CONCLUSIONS 

We present a self-consistent model for simultaneously fit¬ 
ting time-averaged and lag-frequency spectra in the high 
flux state of Mrk 335. To reduce the number of free parame¬ 
ters in our initial model, we restrict our consideration to the 
lamp-post geometry with a maximally rotating black hole 
with a standard disc extending between r ms and 400r s . The 
model including full-dilution effects is successfully fitted to 
the combined 2-10 keV spectrum and the time lags of 2.5-4 
vs. 4-6.5 keV bands. We find the inclination angle is 45°, the 
iron abundance is 0.5 and the photon index of direct con¬ 
tinuum is 2.4. The X-ray source height is very small, ss 2r g , 
illuminating a constant density disc (p = 0). The ioniza¬ 
tion parameter is lO' 5 erg cm s _1 at the innermost part of 
the disc and decreases further out. The black hole mass is 
« 1.3 x 10'M 0 . Our results also show that, along with the 
best value of p , the model still provides excellent fits but the 
constrained source height and black hole mass are different. 
This implies that the well-define density profile is needed in 
order to robustly determine the source height and the cen¬ 
tral mass. Having successfully fitted a physical reverberation 
model to the spectrum and time lags in the hard band, we 
discuss the difficulties of fitting the soft band due to the soft 
excess since its origin is unclear. Nevertheless, we have pre¬ 
sented a plausible model that successfully fit this band by 
further adding only the components which have no (or less) 
effects on reverberation lags. These components, unblurred 
xillver and vmekal, confirm an existence of both cold and 
ionised gas at very distant. In the Appendix, the mismatch 
of soft-excess between xillver and reflionx models is also 
reported which makes the interpretation of the physics at 
0.3-2 keV much more complicated. This is however worth 
investigating in the future. 

Furthermore the positive lags which we are modelling 
phenomenologically by a power-law could be replaced by the 
hard lags from a physical model. The choice of mechanisms 
that drives positive lags (e.g. propagating fluctuations on 
the disc) is also important. Modelling the propagation of 
accretion rate fluctuations and an associated spatially ex¬ 
tended corona geometry is planned for a future paper, but 
is significantly beyond the scope of this paper. This is par¬ 
ticularly important for constraining the otherwise fairly ar¬ 
bitrary power law that can be applied to the frequency de¬ 
pendent time lags. 

The model could also be improved by allowing the outer 
radius to be a free parameter as it affects reverberation lags 
(see Fig. [4|. Additionally, we plan to include the energy- 
dependent time lags in our fitting procedure. All are possi¬ 
ble and straightforward but, undoubtedly, make the model 
more computationally intensive. Fitting combined spectral 
and timing data is important as the parameters are better 
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constrained than fitting each one alone. It also has the po¬ 
tential to constrain or measure the spectrum of reflection 
models in the soft X-ray band. 

Finally, there are two spectral components that our 
model does not capture - reflection from distant neutral 
material, and thermal emission from some hot gas whose ori¬ 
gin is less clear. Both additional components are needed in 
“traditional” models of the broad iron line in time-averaged 
spectra. Future observations with micro-calorimeters will re¬ 
ally help constrain the location and dynamics of the gas as¬ 
sociated with these narrow lines. The features from distant 
reflection, however, do not vary on the timescales of the in¬ 
ner disc reflection which make the reverberation technique 
very powerful. 

We have demonstrated the feasibility of successfully fit¬ 
ting cross spectra to constrain the corona and disc geometry. 
These initial fits reveal some degeneracies between param¬ 
eters that can be investigated by looking at different flux 
states and by analysing a sample of objects. A more sophis¬ 
ticated, and likely spatially extended, model must be de¬ 
veloped to self-consistently explain the positive lags at low 
frequency. For example, at the moment it is far from clear 
how the propagating fluctuations connect to the corona. The 
“soft excess” continues to warrant further investigation. 
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Figure Al. Residuals in the 0.3-2.0 keV band for the REVB ® 
REFLIONX (top panel) and REVB (g> xillver (bottom panel). Each 
model has been fitted from 2—10 keV adding only two narrow 
Gaussians at Ri 6.4 keV (Fe Ka) and 7 keV (Fe XXVI) and then 
extrapolated to lower energies. 



P 

Figure A2. Best fit values of source height (top panel) and black 
hole mass (bottom panel) for variable disc-density index, p (n oc 
r~P). 


parameters. The black hole mass is « 7.9 x 10 6 Mq for the 
reflionx model and « 2.0 x 10 7 Mq for the xillver model. 
The source’s height using reflionx model is 7 r g which is 
much larger than in Emmanoulopoulos et al. (2014). On the 
other hand, the xillver model places the source much lower 
at 2r 9 . Changing the disc density profile results in a different 
disc ionization profile, which in turn produces very different 
reflection spectra and time lags for the two reflection models. 
The disagreement between the two models, especially when 
p = 0, confirms that in some cases the source height and 
the central mass are model-dependent. However, only the 
case of h — 2 r g , as shown in Section [4] can provide the best 
simultaneously fits the time-averaged spectrum and the lags 
in both soft and hard bands. 


APPENDIX 


For comparison, and to assess the systematic differences be¬ 
tween different reflection models, we couple the revb model 
to the reflionx model using similar grid cells as shown 
in Table [T] The revb 0 reflionx provides excellent fits 
to the 2-10 keV spectrum and the hard band la gs w ith 
X 2 /d.o.f. = 0.948. However, the soft excess (Fig. Al I is 
significantly different to that of the revb 0 xillver. The 
parameters which have successfully fitted the soft excess 
for xillver model (Table [3| then should not fit for the 
reflionx model. The mismatch of soft-excess can be under¬ 
stood in term of differences in their atomic data if the xil¬ 


lver model is more absorbed at low energies (Garcia et al. 


2013). These systematic differences are certainly larger than 


our statistical error bars and can affect the interpretation of 
spectra and time lags. 

Keeping in mind the uncertainty of physics below 2 keV, 
we only compare both models by considering the 2-10 keV 
spectrum and the hard band lags. It is interesting to note 
that the best-fit values of the source height and the black 
hole mass significantly depend not only on the disc density 
profiles, but also on the choices of reflection models (see 
Fig. A2). Within the grid cells we investigate, a constant 


density disc (p = 0) had the largest difference in model 
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